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Abstract. Our in silico model was built to investigate the development pro¬ 
cess of the adaptive immune system. For simplicity, we concentrated on hu¬ 
moral immunity and its major components: T cells, B cells, antibodies, in¬ 
terleukins, non-immune self cells, and foreign antigens. Our model is a mi¬ 
croscopic one, similar to the interacting particle models of statistical physics. 
Events are considered random and modelled by a continuous time, finite state 
Markov process, that is, they are controlled by independent exponential clocks. 
Our main purpose was to compare different theoretical models of the adaptive 
immune system and self-nonself discrimination: the ones that are described 
by well-known textbooks, and a novel one developed by our research group. 
Our theoretical model emphasizes the hypothesis that the immune system of a 
fetus can primarily learn what self is but unable to prepare itself for the huge, 
unknown variety of nonself. 

The simulation begins after conception, by developing the immune system 
from scratch and learning the set of self antigens. The simulation ends several 
months after births when a more-or-less stationary state of the immune system 
has been established. We investigate how the immune system can recognize 
and fight against a primary infection. We also investigate that under what 
conditions can an immune memory be created that results in a more effective 
immune response to a repeated infection. 

The MiStlmm simulation software package and the simulation results are 
available at the address http://kerepesi.web.elte.hu/MiStImm/ 


1. Introduction 

Germain et al. [T. Table 1] broadly classifies the computational models in im¬ 
munology into four groups: individual particle-based stochastic, particle number 
stochastic, concentration-based spatial non-stochastic, concentration-based non- 
spatial non-stochastic (see 330 references therein). Our computational model best 
fits into the second category: particle number stochastic. A great advantage of 
such a model is that it can easily incorporate the most important types of cells 
and molecules together with their essential features and events that play important 
roles in immune reactions. In such a model events - for example interactions of 
components - occur at random. Also, such a model is typically microscopic in space 
and limited to a small variety of cells and molecules. 

To further simplify things, we chose the humoral immune system as the first 
modeling objective, since the humoral phase (in particular: blood) may be consid¬ 
ered spatially homogeneous; thus a microscopic spatial volume may represent the 
whole humoral phase well. A major advantage of this approach is that it is not 
necessary to describe the actual spatial positions and spatial motions in the model. 
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Instead, model components randomly choose one of the other components as inter¬ 
action partners, because any components are close enough to become engaged in 
an interaction. 

One of the major objectives that we aimed at was to compare the more-or- 
less generally accepted ways of operation of the immune system as described in 
standard textbooks and a partially new idea of operation that our group had been 
proposing for some time then mm- For brevity, the more-or-less generally 
accepted theoretical models 1113 will be called Conventional Role of Self models 
(CRS models), while our group’s theoretical model will be called the Enhanced Role 
of Self model (ERS model) from now on. For, our theoretical model emphasizes 
the hypothesis that the immune system of a fetus can primarily learn what self is 
but unable to prepare itself for the huge, unknown variety of nonself. 

2. The ERS theoretical model 

2.1. Arguments that led to a novel immune model. 

2.1.1. How does the immune system works? Recently, Thompson j8] painted a pic¬ 
ture of the human immune system with broad brush strokes, which represents the 
consensus view that the immune system is a complex and powerful defense mecha¬ 
nism. We have a somewhat modified picture, Figure |T] 

Tauber, in contrast, suggested that “... ‘immunity’ may be a semantic trap that 
has confined our understanding of the immune system to only a narrow segment of 
defensive, aggressive functions” [5]. 

Trying to resolve this issue, we suggested a paradigm shift in the immune surveil¬ 
lance theory postulating that the evolutionary pressure driving the creation of the 
T cell receptor (TCR) repertoire was primarily the homeostatic surveillance of the 
genome [I]. This is achieved by positively selected T cells, which form a homeostatic 
coupled system via internal dialogue with tissue cells. An important contribution 
was the complementarity theorem of Dillon and Root-Bernstein uni mi, which 
states that a dynamic steady state is achieved through low affinity complementary 
TCR-MHC interactions between T cells and host cells m3- Homeostasis, which 
includes immunity, is controlled by a so-called coupled system via internal dialogue 
between tissues and T cells, while molecular complementarity puts strict limits on 
variations. Therefore, in our view, the primary function of the immune surveil¬ 
lance is homeostatic such that T cells prevent the disappearance of differentiation, 
the natural tendency of people to turn into tumors. Such shift of emphasis from 
infections to cancer seems to be consistent with epidemiologic observations [5]. 

Based on information theoretical principles and the law of parsimony we sug¬ 
gested earlier that recognition of self peptides is sufficient to attack nonsclf mm 
To this end, we proposed that the positively selected subtly variable T cells consti¬ 
tute an evolutionary link between the invariable innate and the hypervariable B cell 
systems. In order to discriminate self and nonself, T lymphocytes should primarily 
recognize the much smaller and always available set of self antigens, rather than 
the practically unlimited and for the immune system only partially known nonself 
antigen universe. 

Thus far, our theoretical models primarily described the role of T cells. The 
objective of the present work is to supplement it with the affinity maturation of B 
lymphocytes in which random events are perhaps the most characteristic. 
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Figure 1. How immunity works: Innate immune cells engulf and 
kill pathogens and release cytokines to enhance the immune re¬ 
sponse. Cytotoxic T cells kill pathogens and attack infected cells. 
B cells and helper T cells interact via interleukins and directly 
to instigate an enhanced immune response. Some hypermutated 
B cells become plasma cells and make antibodies that target spe¬ 
cific antigens. Some B and T cells become memory cells that can 
quickly fight future infections by the same pathogen. 


2.1.2. The homunculus concept. For better understanding our B cell model, a short 
recollection of the homunculus concept and the impetus for the development of our 
one-signal T cell model will be useful. 

To our best knowledge the immunological homunculus concept of autoimmunity 
theory, which was developed by Irun Cohen mmmm, is the closest to our 
current model. It states that the immune system is composed of networks of inter¬ 
acting cells and molecules the aim of which is not to discriminate self from nonself, 
but continuously respond to self in order to organize inflammation in a way that 
maintains, heals and regenerates damaged cells and tissues. This is carried out with 
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a high frequency of lymphocytes that recognize key self-antigens: the immunolog¬ 
ical homunculus. The immunological homunculus positively selects T cells and B 
cells responding to key body molecules forming a functional “internal image” of 
the body. The evolution of the immune system has provided a multilevel system 
that interconnects the innate and adaptive immune systems to serve at least three 
central purposes: the defense from microbial pathogens, the capacity for discrimi¬ 
nation of self from nonself necessary for the prevention of autoimmune disease, and 
essential effector roles in wound repair and tissue remodeling. The immune system, 
therefore, is a cognitive system. 

2.2. Main features of the ERS theoretical model. 

2.2.1. A single T cell cannot discriminate self and nonself, only a wide T cell reper¬ 
toire can. One of the most important points in our argument is that a single T cell 
alone is unable to distinguish between self and non-self, only a wide T cell repertoire 
can. 

Shapes of self and nonself entities are intricately interwoven sets; in the lan¬ 
guage of the shape space model (see below), the subsets of points representing self 
and non-self are complexly interlaced and cannot be separated by a nice smooth 
mathematical curve. Therefore the complexity of the antigen universe exceeds the 
capacity of an individual T cell. The “knowledge” of each specific T cell is reflected 
by the shape of its TCR. An individual T cell therefore is able to recognize only 
a complementary or near complementary MHC-peptide molecule. This way, any 
given T cell clone is responsible for a tiny set of nearly complementarily shaped self- 
MHC-peptide complexes. In the present paper T cells with nearly complementary 
TCR to self-MHC-peptide complexes are designated as regulatory T cells. 

In particular, the complete repertoire of regulatory T cells (see in particular 
Fig. [3] below), is able to reflect the whole set of self antigens. The repertoire of 
regulatory T cells is first created in the thymus of the fetus by negative and positive 
selection and it constitutes the basis for self-nonself discrimination. Any self-MHC- 
peptide complex that is able to attach to one of the regulatory T cell repertoire with 
intermediate affinity can be classified as self ; any other MHC-peptide complexes 
- that has weak affinity to each regulatory T cell or strong affinity to one of the 
T cells - can be classified as nonself. Thus the regulatory T cell repertoire - like 
the conductor of an orchestra - controls other elements of the adaptive immune 
system. This does not exclude the possibility that regulatory T cells - like players 
of an orchestra - may take part in immune reactions similarly to other T cells as 
well. See further details in |5J. After birth, development of infection specific T cell 
and B cell clones are under regulatory T cell control, see below. 

2.2.2. Different T cell - B cell interactions in the serum. As in our current in silico 
model T cell - B cell interactions in the humoral phase is a basic phenomenon, here 
we describe three different types of it. 

In a healthy individual during intrauterine life, randomly produced moderately 
self-reactive B cell clones are confronted with an overwhelming quantity of self 
antigens. Those B cells (and APCs) that attach with intermediate affinity to any 
of these self antigens will present self peptides in their surface MHCII molecules to 
regulatory Th cells. This ensures B cell and regulatory Th cell survival, respectively, 
but it is insufficient to trigger extensive clonally based B cell expansion required for 
specific immunity or autoimmunity. It will be called weak interaction and division 
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from now on. Thus positively selected regulatory Th cells are critical parts of 
the homeostatic control in our model, so that regulatory Th cell clones exist for 
practically all kinds of self-MHCII - self-peptide complexes presented by any of 
the APC’s. After birth, this process maintains an immune image of self which can 
control self-nonself discrimination. 

During an infection a new antigen (e.g. bacterium) appears in the serum. Some 
APC’s, in particular B cells, with appropriate affinity for the new antigens, present 
foreign peptides on their surface MHCII proteins. Since in our model foreign pep¬ 
tides transiently inhibit the complementary TCR-MHC interactions in general, such 
perturbation creates steric hindrance that “obstructs” the docking of positively se¬ 
lected regulatory Th cells to their cognate ligands. Disruption of such contact 
between B cells and regulatory Th cells for a critical period of time by a foreign 
peptide results in an emergency. In order to reestablish contact, foreign peptide 
presenting B cells will secrete chemotactic alarm signals (“smoking gun”) attract¬ 
ing Th cells to this region. This initiates a non-specific, polyclonal activation in 
local Th lymphocytes via the CD28 receptor alone )16| such that a limited beneficial 
local cytokine storm is generated in local Th cells triggering B cells to clonal expan¬ 
sion, hypermutation, and eventually they develop into specific antibody producing 
plasma cells. This will be called intermediate interaction and division from now on. 
The resulting inner state of the affected Th and B cells will be called “activated” 
state. 

The default mode of our model is that a random peptide decreases complemen¬ 
tarity between a naive TCR and the MHC. However, following the initial polyclonal 
activation phase, there is always a possibility that rare T cell and B cell clones with 
higher affinity may well recognize foreign antigens, particularly when a significant 
fraction of host cells is infected and viral load is high (for example in hepatitis, 
see in El). Such higher affinity interactions would then drive clonal (e.g. HCV 
specific) T cell proliferation, activation, lysis of infected cells, as described by the 
conventional two-signal models. Having cleared the infection, specific T cells could 
eventually become an expanded memory type T cell clone, while B cells could dif¬ 
ferentiate into antibody producing plasma cells. It is thought that acquisition of 
memory T-cell function is an irreversible differentiation event. Unlike regulatory 
T cells, such population does not require self-peptide-MHC complexes for mainte¬ 
nance. Nevertheless, sustaining the functional phenotype of T memory cells requires 
active signaling via CD27 m- Specific T and B cell activation, proliferation and 
lysis of infected cells, therefore, obey the rules of the conventional two-signal model. 
Clearly, this process may require several days in general. It will be called strong 
interaction and division in the sequel. The resulting inner state of the affected Th 
and B cells will be called “strongly activated” state. 


3. Description of the MiStImm in silico model 

3.1. Basics. We call our in silico model Microscopic Stochastic Immune model or 
briefly MiStImm model. It is a further developed version of our 1994-98 B cell 
model [19]. According to the broad classification by Germain et al. [IJ Table 1] 
our model is a particle number stochastic computational model. Our software is a 
C program, though it was written in the spirit of the agent-based models. 
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Theoretically, the interactions of the components and other events are described 
by a continuous time, finite state Markov process. A Markov process is a memory¬ 
less stochastic process: if we specify the present state of the system, then we may 
forget about its history when we want to investigate its behavior in the future. 

Our model has finitely many components: helper T cells (regulatory Th cells and 
potentially infection/mutation specific Th cells), B cells, antibodies, interleukins, 
non-immune self cells, and foreign antigens. Presently, other than helper type T 
cells or other antigen presenting cells besides B cells are not represented in our com¬ 
puter model. Each component has a number of characteristics (parameters) and 
certain attached random events or processes of events that may occur at random. 
This means that to each potential event there belongs an independent random expo¬ 
nential clock. That is, the model is simulated by independent exponential random 
variables. The expected waiting time r of each clock may depend on the present 
state of the whole system. The main step of the algorithm in the computer model 
is that the event whose clock rings first occurs. For example, the corresponding 
component interacts with a randomly chosen partner. The occurrence of such an 
event may cause several changes in the model, like births, deaths, and updates of 
parameters. Also, in general, expected waiting times of several clocks have to be 
updated. Then the Markov process - that is, each existing clock starts again with 
the new settings. 


3.1.1. Some related computational models. One of the first experiments with a de¬ 
tailed agent-based model (IMMSIM) of immune system was the work of Celada 
and Seiden nn ns [22]. Their goal was to capture the dynamics of the immune 
system and to perform experiments in silico. In [23], they studied the thymus, the 
regulation of positive and negative selection, and the dynamics of the production 
of the TCR repertoire in the thymus. 

A closely related agent-based model, the C-ImmSim package has been developed 
and investigated by Rapin et al. An excellent recent description of their work can be 
found in |24| . Their model represents pathogens, as well as lymphocytes receptors, 
by means of their amino acid sequences and makes use of bioinformatics methods 
for T and B cell epitope prediction. This is a key step for their simulation of the 
immune response, because it determines immunogenicity. The related book 125] 
can be used as a practical guide to implement a computational model with which 
one can study a specific disease. 

The Basic Immune Simulator (BIS) [26] is also an agent-based computing model 
to study the interactions between innate and adaptive immunity. The BIS was 
created using the Recursive Porus Agent Simulation Toolkit (RepastJ) library, an 
open-source software library that is available online. 

Finally, Kalita et al. m have developed SIMISYS, which is also a cellular 
automata model of the human immune system. It uses tens of thousands of cells 
and innate and adaptive components of the immune system. In particular, the 
model contains macrophages, dendritic cells, neutrophils, natural killer cells, B 
cells, T helper cells, complement proteins, and pathogenic bacteria. 


3.1.2. Peptide lattice. Our computer model takes a microscopic volume of the hu¬ 
moral phase and also a microscopically small part of the shape space universe. To 
explain what we mean by shape space, assume that the shape of a T cell receptor 
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Figure 2. Two simplified complementary shapes characterized by 
the points ( Xp,yp ) and (xp, —yp), respectively 


(TCR) can be represented by a point in a large set of a Euclidean space. Theo¬ 
retical considerations compared with experimental data led to the conclusion | 28 j 
that the dimension of this shape space, i.e. the number of parameters essential in 
describing a binding, is not too large, probably around five. 

The microscopically small part of the shape space that we consider in our model 
is a small discrete N x N planar grid in the shape space (e.g. N = 1000). The x £ 
{0,1,..., N} coordinate of a shape point may represent a “horizontal” coordinate of 
the main part of the binding profile of a TCR or an MHC+peptide complex, while 
the y £ {— N/2 ,..., N/2} coordinate may represent the “vertical” coordinate of 
the main part of the binding profile. A positive coordinate represents “convexity”, 
while a negative coordinate represents “concavity”. Figure [2] shows our underlying 
idea for the shape of a peptide characterized by a single point ( xp,yp ). Needless 
to say that our model of shapes is a much simplified one, but is still suitable to 
represent essential binding properties of antigens. We call the above finite square 
grid the peptide lattice in the sequel. 


3.1.3. Antigen lattice. Shape of a B cell receptor (BCR) or shape of an antigen 
is similarly represented by a point of an antigen lattice in the model. Here again 
the x £ {0,1,...,TV} coordinate of a shape point may represent a “horizontal” 
coordinate of the main part of the binding profile of the BCR or antigen, while the 
y £ {—N/2,...,N/2} coordinate may represent the “vertical” coordinate of the 
main part of the binding profile; a positive coordinate representing “convexity”, 
while a negative coordinate representing “concavity”, see Fig. 1. 

For simplicity, to each antigen (xa,Va) in the antigen lattice we assign exactly 
one peptide (xp,yp) in the peptide lattice. To make identification of an antigen 
and its corresponding peptide easier, we will use the convention that Xa—Xp and 
Va = yp- 

3.1.4. Complementarity. Complementarity plays a basic role in binding. The per¬ 
fect fit between a TCR and an MHC+peptide complex means in the model that 
the shape (xp, yr) of the TCR and the shape (xp, yp) of the MHC+peptide satisfy 
the equalities Xp = xp and yp = —yp, see Fig. 1. In the model we introduce a 
metric or distance d to measure the degree of similarity of two shapes z\ \= (a?i, S/i) 
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Figure 3. Simplified graphical representation of the difference be¬ 
tween the ERS and the CRS model 


and 2 2 := (^ 2 , 1 / 2 ): 

d(z 1 ,z 2 ) := max -{>2 - aq|, \y 2 - Vi\}- 

A TCR zt := ( Xp,yr) and an MHC+peptide Zp := ( Xp,yp ) is nearly comple¬ 
mentary in our model if the distance of zp and Ip := ( xp , — yp) is small enough. 
Similar is the representation of the complementarity between BCRs and antigens 
in the model. 

Only complementary or nearly complementary shaped ligands and receptors can 
bind. The dots in Figure [3] represent TCRs that are exactly complementary to 
some self MHC+self-peptide complex. The areas shaded in darker green are called 
the characteristic rings of self-peptides. They represent the set of shapes that 
are allocated to possible regulatory T cells after negative and positive selection in 
the ERS model, see |3.3.7| below. The areas denoted by lighter green correspond to 
possible shapes of classical, potentially infection (or mutation) specific T cells , while 
white areas are representing self-reactive T cells that are prohibited for T cells in 
the two respective models. Observe that in the ERS model, moderately self-reactive 
T cells are present after negative and positive selection. In fact, they constitute 
the most important class of T cells that decide self-nonself discrimination. On the 
other hand, such moderately self-reactive T cells are negatively selected out in CRS 
models. 

We mention that with the above metric, “circles” are in fact squares in our shape 
space model. 


3.1.5. A logistic function. In biology it is typical that when the size of a certain cell 
population gets larger the per capita birth rate in the population decreases. Thus 
the size of a population first increases fast, later it slows down, and in many cases 
at the end it gets relatively stable, depending on model parameters. So to control 
the birth rate or other quantities we use a class of logistic functions, previously 
applied by many other authors (see e.g. [25JSBD- I n general, the next function is 
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32,3 0*0 32,6 0*0 




Figure 4. Two examples of a logistic function 

a useful tool to smoothly and flexibly bound a quantity: 

(3.1) g e , v (x) := p, + xV = ( X + (j) ) (* > 0; 0 > 0, g > 0). 

This formula describes a decreasing function which is equal to 1 for x = 0, 1/2 for 
the threshold value 9, and goes to 0 as x —1 00 , see Figure |4j Its parameters 9 and 
g are set from case to case. 


3.2. Self cells. At time zero, there is a number (say, 3) of different types of non- 
imrnune self cells (briefly: self cells), each with a given initial population size (e.g. 
150). A certain type of self cells is represented by its position ( xs,ys ) in the 
antigen lattice and its peptide ( xp,yp ) in the peptide lattice. Specifically, there is 
a population of bone marrow cells, handled separately from other self cells, with a 
given initial population size. 

Each type of self cells comes with a birth process with a given initial rate (that 
is, with a given initial average waiting time r s 0 between divisions). If the size of the 
population of a specific self cell at a certain time t is s = s(f), then the conditional 
expected waiting between two divisions in this population is 


(3.2) 


T s = 




sge, v {s) 


= -( 1 + G) 


Formula (3.2) indicates that when the number s of a type of self cells becomes 
significantly larger than its threshold value 9 its division rate gets close to zero. 

The case of bone marrow cells is special because it comes not only with a birth 
rate, but, with given rates, bone marrow cells also produce naive B cells and Th 
cells. Naive B and Th cells have randomly determined BCR and TCR shapes that 
are uniformly distributed on the antigen and peptide lattices, respectively. 


3.3. Th cells. While different types of non-immune self cells and foreign antigens 
are treated as populations, B and Th cells are handled individually in the model. 
Naive Th cells are born in the bone marrow. The birth of a Th cell initiates its 
own (natural) death event and a Th cell action process. 

3.3.1. Th cell recognition region. Each Th cell has a recognition region in the pep¬ 
tide lattice. If a TCR is described by the point ( Xt,Vt ), then the corresponding 
recognition region is a square with center (xt,—Vt) and radius rp- The radius of 
the TCR is a constant, there is no hypermutation or affinity maturation for Th 
cells. The recognition region describes the potential shapes of antigens with which 
a TCR can bind: the smaller the distance between a peptide ( Xp, yp ) located on an 
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MHCII and the center (xp, —Vt ) of the recognition region of the TCR, the better 
the fit. 

3.3.2. Thymus. To each naive Th cell there is assigned a random event that places 
it into the thymus. Here the Th cell goes under a negative and a positive selection 
process. Negative selection kills Th cells that are closer to one of the self-peptides 
than a minimum radius r m i n ; negative selection occurs with a given large proba¬ 
bility, tipically pn = 0.99. 

Positive selection kills Th cells that are farther from each self peptide than a 
maximum radius r max ; positive selection occurs with a given, relatively smaller 
probability, typically pp = 0.9. This way, some of the randomly generated Th cells 
that cannot bound self-pepticles may still survive and they can become infection or 
mutation specific Th cells later. 

The degree of maturity of a naive Th cell is 0. If a TCR is in the characteristic 
ring around the reflected image of some self-peptide (see Fig.2), that is, r m j„ < 
d(zp,zr) < r max , where zp := (xp,—yp) is the center of the recognition region 
of the Th cell and zp := ( Xp,yp ) represents the shape of a self-peptide, then it is 
called a regulatory Th cell and its degree of maturity is set to 2. In our model a 
regulatory Th cell has double role. On one hand, it takes part in the controlling 
role of the regulatory T cell repertoire, but it can also act as a Th cell. 

Other Th cells that have survived the negative and positive selections, but are 
outside of the characteristic ring of each self-peptide, are called potential infection 
or mutation specific Th cells, and their degree of maturity is set to 1. 

3.3.3. Th cell actions. For each Th cell, there is a sequence of actions, with expo¬ 
nential random waiting times between two actions. At each action the Th cell is 
to randomly choose one of the potential target MHCII+peptide complexes in its 
recognition region. The closer an MHCII+peptide complex to the center of the 
recognition region, the bigger its chance of being selected. 

3.3.4. Th cell activation control process. It is a sequence of frequently occurring 
random events whose purpose is to check and possibly change the state of activation 
of a Th cell. A Th cell can be in a state of activated or non-activated. This process 
checks if this Th cell has received interleukin of type 1 in a critical period of time 
before this check. If the result of this check is “y es ”> then the Th cell is set to 
“activated” (stress=l); otherwise it is set to “non-activated” (stress=0). 

An “activated” Th cell starts an interleukin of type 2 sending process. This 
process is a signal of its activated state for “activated” B cells in its environment. 
We use the symbolic names “interleukin of type 1 or 2” in this paper, without 
specifying the exact type of these interleukins. 

An “activated” Th cell begins cell division of intermediate kind. Division of the 
intermediate kind is different from the weak or strong kind, see below. 

3.3.5. Self-nonself discrimination. When a regulatory Th cell (degree of maturity 
is 2) bounds with intermediate affinity an MHCII+peptide complex which has state 
“non-activated” (see below), then with high confidence it means that the peptide is 
a self-peptide. This contact initiates a division of weak kind for both this regulatory 
Th cell and the attached B cell. This weak division helps to stabilize this interac¬ 
tion among three partners: self-cells, B cells that can react to self, and regulatory 
Th cells that can attach to this self-peptide with intermediate affinity. It is also 
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important that B cells that can process self antigens cannot start a intermediate or 
strong division process. This is an important inhibitory effect of regulatory Th cells. 
For, such B cells that react to self are in a state of “non-activated” permanently 
with large probability. 

When a Th cell that has already went through the thymus, but it is not a 
regulatory Th cell (thus its degree of maturity is 1), obtains interleukin of type 1 
then it may begin a non-specific division of intermediate kind and may start to 
secrete interleukin of type 2 to start division of intermediate kind of activated B 
cells. 

If a Th cell has already went through the thymus, but it is not a regulatory 
Th cell (thus its degree of maturity is 1), the target is an activated B cell, and 
the distance of attachment satisfies d(zp,zp) < then with high confidence it 
means that the peptide is foreign or mutated self. Here Jp = (xt, ~Ut) is the center 
of the recognition region of the Th cell, zp := (xp, yp) is the point representing the 
peptide, and r m i n is the inner radius of the characteristic ring around the reflected 
image of self-peptides. Remember that because of the negative selection, such short 
distance between a self-peptide and the center of recognition region is extremely 
unlikely. Then both this B cell and Th cell are very likely useful tools to fight 
against an infection. As a result, this interaction may initiate a division of strong 
kind both in the affected B and Th cells, plus stimulates the secretion of interleukin 
of type 1 (in the B cell) and type 2 (in the Th cell). Strong division of a B cell 
implies its hypermutation with given probability as well. This is a direct help of 
the Th cell for the affected B cell. 

3.3.6. Th cell divisions. The probability of division of a Th cell may depend on 
several factors. It may get bigger when the distance d between the MHCII+peptide 
complex and the TCR is smaller (i.e., the complementarity is better). It gets smaller 
when the number no of all TCR’s is large (i.e., the concentration of Th cells is 
already large). It gets smaller when the number m of TCR’s in a neighborhood of 
the Th cell is large (i.e., the local concentration of Th cells is already large). The 
formula for the probability of division is given by somewhat different formulas for 
the three different kinds of division: weak, intermediate, and strong; namely, not 
each of these depends on all the three factors. 

The probability of a division of weak kind of a Th cell is given by 

(3-3) VT.W — k w 9 $ n 0 , ri n o (^o) QOnhVnl (^r)‘ 

The purpose of division of weak kind is to establish a stable contact between self 
antigens, B cells reacting to self with a weak affinity, and Th cells reacting to self 
peptides with a intermediate, standard affinity. 

The probability of a division of intermediate kind of a Th cell is given by 

(3-4) p T>m = k m ge n0tVn0 (n 0 ) gg nl l7Jnl (rii)• 

The purpose of division of intermediate kind is to create a fast, non-specific im¬ 
mune reaction to a new, typically quickly growing number of nonself antigens. The 
growing amount of Th cell help (interleukin of type 2) can help the division of 
intermediate kind of B cells that are able to bind the new nonself antigens in the 
humoral phase. 

The probability of a division of strong kind of a Th cell is given by 

( 3 - 5 ) p T ,s = K ge dtVd (d) ge n0 ,rino( n o) ge n uvm( n i)- 
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The purpose of division of strong kind of Th cells is to initiate a strong immune 
reaction when infection or mutation specific Th cells appear and can bind infection 
or mutation specific B cells. Important requirements to such a division that the 
binding distance satisfy d < and the attached MHCII be “activated”. These 
requirements can guarantee with large probability that this strong reaction is not 
arising against self. Then this Th cell becomes “strongly activated” (stress=2). 
This condition is independent of interleukins of type 1. 

3.3.7. Regulatory Th cells. As we saw above, the regulatory Th cell repertoire plays 
a most important role in our model. This role is similar to the one they have in 
the computational model [311 . Starting in the fetus, and throughout the entire life 
span, they give a faithful mirror-image of the self-peptide repertoire. 

• They regularly visit B cells having only self-peptides on their MHCII and 
inhibit their strong division, but support their weak division. 

• They are players in normal Th cell roles, like helping non-specific interme¬ 
diate type and specific strong type division of B cells. They can also secrete 
interleukins of type 2. 

3.4. Interleukins. We use the symbolic names “interleukin of type 1 or 2” in this 
paper, without specifying the exact type of these interleukins, similarly to Figure 

3 of [sg. 

Interleukins of type 1 are emitted by activated B lymphocytes. This process 
initiates an action process and also a death process of these interleukins. Each 
interleukin of type 1 molecule randomly chooses a Th cell. This is a signal for the 
Th cell to start intermediate type division and to secrete interleukin of type 2. 

Interleukins of type 2 are emitted by Th lymphocites. This process initiates an 
action process and also a death process of these interleukins. Each interleukin of 
type 2 molecule randomly chooses a B cell. This is a signal (an indirect help) for a 
activated B cell to start cell division of intermediate kind. 

3.5. B cells. Naive B cells are born in the bone marrow. The birth of a B cell 
initiates its own (natural) death event, B cell action process, and B cell activation 
control process, each with separate rate. Each B cell carries a number (say, 3) of 
MHCII molecules. 

3.5.1. B cell recognition region. Each B cell has a recognition region in the anti¬ 
gen lattice. If a BCR is described by the point ( Xb,Vb ), then the corresponding 
recognition region is a square with center (xb,—Vb) and radius rs- The radius of 
the BCR of a naive B cell is a given constant, while B cells that are born in the 
periphery after hypermutation may have smaller radii. The BCR z' B = ( x' B ,y' B ) of 
a hypermutated B cell offspring is determined at random, uniformly on a square 
around the mother BCR. Thus there is only a chance that its affinity to a given 
antigen za = {xaiVa) Is higher , that is, the distance d(zA,z' B ) is smaller than that 
of its mother cell. The radius r' B of a hypermutated offspring will be smaller than 
that of its mother cell depending on the above distance: r' B = crs r$. Typical 
values are c = 0.9 and ro = 5. This effect may increase the affinity of some “lucky” 
offspring to the given antigen. 

In sum, the recognition region describes the potential shapes of antigens with 
which a BCR can bind: the smaller the distance between an antigen za = ( xa , Va) 
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and the center zb = (ib, —j/b) of the recognition region, the better the fit between 
the antigen and the BCR. 

3.5.2. B cell action process. For each B cell, there is a sequence of actions, with 
independent exponential waiting times between two actions. At each action the B 
cell is to randomly choose one of the potential target antigens in its recognition 
region. A target can be another B cell, an antibody, a non-immune self cell, or a 
foreign antigen. The closer an antigen za = {xa, Da) to the center zb = {xb, —Ub) 
of the recognition region, the bigger its chance of being selected as the next target. 
The chosen target can be killed only if the above distance is smaller than the recog¬ 
nition radius tb of the B cell, that is, d(zA,ZB) < Tb- The smaller this distance, 
the larger the probability that the antigen will really be destroyed. Since smaller 
distance represents stronger affinity in the model, it means longer attachment be¬ 
tween an antigen and the BCR. So this condition is equivalent to the fact that a 
target can be killed if it is bound to the BCR for a long enough time. 

When the chosen target is destroyed, its peptide is placed on one of the MHCII’s 
of the B cell. The MHCII selected is primarily an empty one; when all of the 
MHCII’s are already loaded, then one of them is chosen at random to replace the 
old peptide by the new one. 

3.5.3. B cell negative selection filter in the bone marrow. To each naive (immature) 
B cell there is assigned a random event that places it into a negative selection filter 
in the bone marrow. Negative selection kills B cells that are closer to one of the 
self-antigens than a minimum radius r m i n f,; negative selection occurs with a given 
large probability, typically p^b = 0.99. 

The degree of maturity of a naive B cell is 0. A B cell that has survived the 
negative selection is called a mature B cell, and their degree of maturity is set to 
1. Only B cells with degree of maturity > 1 can function as normal B cells. 

3.5.4. B cell activation control process. It is a sequence of frequently occurring 
events whose purpose is to check and possibly change the state of activation of a B 
cell. The main parameter is the critical time t cr it- Each of the MHCII carried by 
a B cell can be in a state of “ activated” or “non-activated”. An empty MHCII is 
not “activated” by definition. 

• A given non-empty MHCII is set to “non-activated” when the time elapsed 
since the last event effecting this MHCII is less than t cr u. Such an event can 
be a regulatory Tli cell attaching to this MHCII, or placing a new peptide 
on this MHCII. 

• A given MHCII is set to “activated” when the time elapsed since the last 
event effecting this MHCII is greater than or equal to t cr it ■ 

Similarly, a B cell can also be in a state of “ activated” or “non-activated”. 

• When its each MHCII is in the state of “non-activated”, the B cell itself is 
set to state of “non-activated”. 

• When at least one of its MHCII is “activated”, then the B cell is set to 
“activated”. 

An “activated” B cell starts an interleukin of type 1 sending process. This 
process is a signal of its activated state for Tli cells in its environment. 

An “activated” B cell may start a cell division of intermediate kind if it obtains 
help from non-specific Th cells. Help may come as interleukin of type 2 produced 
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by Th cells, that has arrived in a critical period of time before this check. (This 
kind of cell division cannot occur with plasma cells or memory cells.) Division of 
the intermediate kind is different from the weak or strong kind. Here the activation 
(stress) level is 1. 

In the case of cell division of the strong kind , which occurs by the help of infection 
or mutation specific Th cells, the activation (stress) level is 2. 

3.5.5. B cell division and maturity. Each B cell has a degree of maturity. A na'ive, 
immature B cell has degree 0, while B cells that have survived a negative selection 
filter in the bone marrow are mature B cells, having degree of maturity 1 first. 
Mature B cells may encounter antigens at the periphery. A B cell division can be 
the result of an encounter with an antigen which is escorted by a direct or indirect 
(via interleukin) help from a Th cell, see below. At each division of a B cell, one 
of the two offspring inherits all characteristics of the mother cell (let us call it the 
first offspring for explicitness), while the other offspring (let us call it the second 
offspring) may undergo hypermutation with given probability. The first offspring 
inherits the mother’s MHCII-peptide complexes, while the second offspring starts 
with empty MHCII’s. The second offspring after the first division has a degree of 
maturity 2. The result of a hypermutation is a B cell with randomly shaped BCR. 
The possible shapes are uniformly distributed on a square of the antigen lattice, 
with given radius around the mother BCR. 

A second division may lead to two different outcomes with given probabilities: 
the second offspring can be either a memory cell (degree=3) or a plasma cell (de- 
gree=4). A memory cell has the same characteristics as a normal B cell except 
that its average lifespan is significantly longer (e.g 10 days instead of the standard 
3 days). A plasma cell constantly - at random time instants - produces antibodies 
of the type of its own BCR. 

Possibility of division of a B cell arises after contacting an antigen or obtaining 
Th help in the form of interleukin of type 2. The probability of division of a B cell 
depends on several factors. It gets bigger when the distance d between the antigen 
and the BCR is smaller (i.e., the complementarity is better), or when the radius r of 
the recognition region of the BCR is smaller (i.e., the affinity of the B cell is bigger). 
It gets smaller when the number no of other BCRs in a rectangle around the BCR 
is small (i.e., the concentration of B cells is already large). Finally, one or two 
factors can depend on the concentration difference c between the number of targets 
in the recognition region of the B cell and the number of targets in the reflected 
image of the recognition region. If the concentration difference is too small, the B 
cell may get insensitive. If the concentration difference is too large, the B cell may 
get anergic. 

The specific formulas for the probability of division in the respective cases of 
weak, intermediate, and strong B cell divisions are as follows. The probability of a 
division of weak kind of a B cell is given by 
(3.6) 

PB,w = k w gr minb +B d ,ri d (d)(l — 9r minb —6 d ,r} d (d)) g8 r ,r),. ( r ) ( n o) (1 — 9n m ,r) a ( c ))- 

The purpose of division of weak kind is to establish a stable contact between self 
antigens, B cells reacting to self with a weak affinity, and Th cells reacting to self 
peptides with a intermediate, standard affinity. The first, constant factor k w is 
typically 1. The purpose of the second and third factors depending on d is to help 
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those B cells that are at a standard distance from their targets, in the present 
case, self antigens. The last factor, depending on c intends to guarantee that a 
large number of antigens, typical for self antigens, be in the recognition region of 
the weakly dividing B cells. The first parameter n m here is the actual number of 
bone marrow cells, which is a common measure of the size of non-immune self cell 
populations. 

The probability of a division of intermediate kind of a B cell is given by 

(3.7) p B ,m = k m g 0 d , Vd {d)ge r ,r lr (r)ge n ,ri n (no)ge C 2 ,vc( c ) (! - 50ci,>?c( c ))- 

The purpose of division of intermediate kind is to create a fast, non-specific im¬ 
mune reaction to a new, typically quickly growing number of nonself antigens. The 
growing amount of B cells that are able to bind the new nonself antigens in the 
humoral phase even when there exist no infection or mutation specific B or Th 
cells can give an early start to an effective immune reaction. Activated B cells can 
release interleukin of type 1 to initiate a non-specific Th help as well. The value of 
the constant multiplier k m is typically 100. 

The probability of a division of strong kind of a B cell is given by 

(3.8) p B ,s = k s ge d , Vd (d)ge r ,r, r (r)g0 n ,ri n (no)g9 c2 ,Ti c (c) (1 ~ ge cl ,nc( c ))- 

The purpose of division of strong kind of B cells is to initiate a strong immune 
reaction when infection or mutation specific Th cells appear and can bind infection 
or mutation specific B cells. Important requirements to such a division that an 
“activated” Th cell binds an “activated” MHCII of this B cell and the binding 
distance between the reflected image of the TCR and the peptide is smaller than 
rm r f n . These requirements can guarantee with large probability that this strong 
reaction is not arising against self. The value of the constant k s is typically 200. 

3.5.6. B cell affinity maturation and network memory. Like in natural selection, 
there exists neither intelligent control which would direct genetic mutations toward 
better fit, nor memory that would save cells from genetically searching a proved 
wrong “direction”. The major effect which has physiological consequences on a B 
cell is the strength of antigen binding. This is like finding the source of heat in 
a dark room, using a single thermometer, with no direct sensing of direction and 
with no memory. The technique the present model applies is a microscopic analog 
of evolution: hypermutation and selection, with survival of the fittest. Namely, the 
program uses a stochastic search for best fit (or a stochastic learning process): 

• An offspring may be randomly hypermutated, so a random variation is 
created in the affinity to the given antigen. 

• The stronger a B cell can bind an given antigen, the more offspring it can 
produce. 

• When the concentration of the given antigen is decreasing, a competition 
arises among B cells for the antigen, and those having higher affinity would 
win in this selection process. 

An affinity maturation model has to handle the danger of autoimmunity. Even 
if naive Th cells which can strongly bind self peptides are deleted as a result of 
negative selection in the thymus, and also naive B cells which can strongly bind 
self antigens are deleted as a result of negative selection in the bone marrow, still 
there is the danger that autoimmune B cell clones may be produced as a result 
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of hypermutation. In the presented model there is a double defense against this 
danger. 


The absence of T cell help in the case of B cells that react strongly to non- 
immune self antigens inhibits their division. This is an essential difference 
between self and nonself in the model. 

Since nonself antigens which can start somatic hypermutation typically 
appear after birth, when the number of self cells is already very large, 
one can argue that at that time randomly produced self-reactive B cell 
clones are confronted with an overwhelming quantity of self antigens. As a 
result, these B cell clones would become anergic [33j. In the model this is 
simulated in the B cell division process: divisions of a B cell, see (3.71 and 
(3.81, become less frequent when the number of objects in its recognition 
region becomes overwhelmingly large. The reproduction process of B cells 
is fastest when the concentration of the complementary antigens is neither 
too small, nor too large. This is common for both self or nonself antigens 
in the model, so when nonself overgrows an upper threshold, the model 
immune system remains practically defenseless against it as well. 


As a result of the double defense described above, there will be “holes” in the 
adaptive immune system, both in the T cell and B cell populations, around the 
mirror image of non-immune self cells l34l l35| . The negative selection in the model 
is especially important during early ontogenesis when the smaller population of host 
cells is vulnerable to self-reactive immune cells. As the individual reaches adult 
size, the large number of host cells plus the absence of T cell help can alone inhibit 
reproduction and affinity maturation of immune cells. Then negative selection in 
the model (like in reality in the thymus) becomes less essential. 

It is reasonable to expect that after a somatic hypermutation - affinity matura¬ 
tion process the resulted specific B cell clones may survive for a certain period of 
time as a local memory. In the model, expansion of certain B cell clones (e.g. as a 
result of an infection by a foreign antigen), under favorable conditions, stimulates 
the reproduction of secondary B cells which are complementary to the expanded 
primary B cell clones and whose receptors are, therefore, similar to the infecting 
antigens. (Of course, similarity here means a mimicry of a binding partner and not 
similarity at the molecular level.) Thus a mirroring process (“ping-pong”) and a 
local network memory may develop and last for a longer time, even in the absence 
of the stimulating antigen. While this memory lasts, repeated infection of the same 
pathogen is eliminated more efficiently. This network model of immune memory 
essentially conforms to Jerne’s immune network concept. Beside other factors, like 
longer living memory cells or antigen preserving follicular dendritic cells, this could 
be a possible explanation of immune memory. 


3.6. Antibodies. A plasma cell is a special kind of B cells, a result of a B cell 
maturity process, cf. 3.5.5 A plasma cell has neither a B cell action process (cf. 


3.5.2), nor a B cell activation control process (cf. 3.5.4). On the other hand, it has 


an antibody birth and an antibody death process. An antibody has the same shape 
in the antigen lattice as the BCR of its mother plasma cell. 

Antibodies have similar action processes as B cells, but, naturally, when destroy¬ 
ing a target, peptide of the target does not appear on an MHCII. The complement 
sub-system of the immune system is currently not represented in the model, so it is 







A MICROSCOPIC STOCHASTIC MODEL OF THE IMMUNE DEVELOPMENT 


17 


supposed that attached antibodies lead to the destruction of the targeted antigen 
with certain probabilities. 

3.7. Foreign antigens. After birth, different nonself antigens may enter the body, 
perhaps several times (e.g. repeated infections with the same disease). A foreign 
antigen is represented by its position (xf,Vf) in the antigen lattice and its peptide 
(: KfpiVfp) in the peptide lattice. A foreign cell comes with an initial population 
size and a birth process with a given initial rate (that is, with a given initial average 
waiting time r/o between divisions). If the size of the population of a specific 
nonself cell at a certain time t is / = /(f), then the conditional expected waiting 
time between two divisions in this population is 



(3.9) 


4. Simulation results 


4.1. Initial parameters. In this chapter we assign specific parameters to the 
MiStlmm in silico model discussed above and we describe the results obtained 
by the corresponding computer simulations. Due to computational limitations, 
simplifications in the parameter settings were necessary. Nevertheless, we tried to 
make the model so that it can reflect the processes that play an essential role in 
the self-nonself discrimination of the adaptive immune system. 

The simulation program can be initialized by about one hundred input parame¬ 
ters, most of which are used to set various immune system models, first of all the 
above mentioned ERS and CRS models. Once one has selected a model, there are 
parameters with which one can make individual specific settings and start various 
event processes (for example different foreign cell injections). Furthermore, we have 
the opportunity to set different random numbers to run different realizations with 
the same parameter settings, thereby observing the role of randomness. 

The following parameter settings were chosen so that they comply with the ERS 
model and if nothing else said we keep these settings. Firstly, we wanted to show 
that a simulation driven by the ERS model can comply with the basic requirements 
that are expected from an immune system model. Secondly, we used simulation 
results of the ERS and CRS models to compare the two. We hope that our ideas and 
our in silico model may encourage investigations about the disputed problems by in 
vitro and in vivo experiments as well. We would especially like to see experiments 
clarifying questions about the self-nonself discrimination. 

4.2. Development and homeostasis of the immune system. The simulation 
starts from a few days after conception and goes until the 5000th time instant, the 
unit of time being a tenth of a day (2 hours and 24 minutes). Initially there are 
only three types of non-immune self cell populations in the model, each with 150 
cells, and no other components. Each of these populations is accompanied by a cell 
division process which results a continuous growing of the number of cells in the 
population, with decreasing rate in time, as one might expect (Figure [5]). 


Coordinates of different antigens of the three different self cells are respectively 
(550,300), (700,-200), (850,150) in the {0,1000} x {-500,500} antigen lattice. 
Coordinates of peptides belonging to these self cells are chosen the same values in 
the {0,1000} x {—500, 500} peptide lattice. 
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Figure 5. Development and homeostasis of immune system 


Bone marrow cells, and the B and T cells generated by them, first appear at 
the 100th time instant. The number of these cells also grow continuously with a 
decreasing rate. TCR rings around the mirror images of self peptides (letters ‘s’) 
- that are characteristic features of the ERS model - begin to develop about the 
1500th time instant and more or less stabilize by the 2800th time instant (Figure|6|. 
There are random fluctuations of these rings: Sometimes the local and global T cell 
populations overgrowth the upper limits and this reduces the probability of T cell 
division. Another reason is that there is a fluctuation in the peptide presentation 
of B cells. 

In sum, one can observe that the qualitative behavior of the development and 
homeostasis in the immune system are not significantly affected by random factors. 

4.3. Normal immune response. In this subsection we show that MiStlmm can 
model normal immune response. Firstly, immune system with normal parameters 
do not attack self cells strongly, only to a very limited extent. Some B cells must 
continuously present self peptides to ensure that T cell characteristic rings around 
self peptides are persistently maintained. Because of anergy, this type of immune 
response is weak and always settles down quickly before it becomes pathological. 

Secondly, a normal immune response should have the ability to destroy the ma¬ 
jority of pathogens - some of them quickly, others perhaps slowly, while in some 
cases it may fail. In our model, death of an individual occurs when the pathogen 
population grows up irreversibly so that its size reaches 4000 items. Diversity of 
foreign intruders are represented by different locations of their receptors and dif¬ 
ferent speeds of growth. In the current settings we work with one kind of antigen. 
Table [l] shows the frequency of winning against pathogen destroyed by a healthy 
immune system - considering three different growth rate and three different initial 
sizes. Figure [7] shows a typical B cell response to an infection over the antigen 
lattice; the position of nonself population is denoted by a letter ‘n’. 

The third important feature of a normal immune response is immune memory. 
Thanks to memory cells, a second immune response against the same infection will 
be more effective than at the primary infection. We have tested this phenomena 
with 100 simulations, adding the same pathogen at 3000th and 3150th time instant. 
There was 81 cases when immune system won against both infections and these 
times the average time interval of the wins were 88.91 and 53.28 (the first and 
the second occasions, respectively). The p -value of the two samples was less then 
0.0001. Figures [8] and [9] show an example for an infection that was repeated once. 
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Figure 6 . Peptide lattice with rings around the reflected images 
of non-immune self cells 


Experiments with our microscopic computational model showed that it cannot 
fight effectively against more than a couple simultaneous infections. Similar is the 
case with the development of immune memory. This observation fits experimental 
results [36]. 

4.4. Autoimmunity. Lack of negative selection of B cells results autoimmunity 
(Figure [l0|. 


4.5. ERS model vs CRS model. One can switch the ERS model to a CRS 
model by modifying three parameters. Simulation results with a CRS model was 
very similar to simulation results obtained with the ERS model above, except the 
establishment of characteristic T cell rings around self antigens. One can say that 
CRS models also satisfy the basic requirements that are expected from an immune 
system model. 

However, the two model are different in one important aspect. We compared the 
efficiencies of the immune reactions in the two models. Our results showed that in 
the ERS model the adaptive immune reaction started faster and was able to destroy 
infections with critically large initial sizes more often than a CRS model (Table [l]). 
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Figure 7. Antigen lattice, B cell response to a pathogen 


Table 1. Normal immune response against of the ERS model vs 
the CRS model simulated by MiStlmm 500-500 times at various 
type of infections. Foreigns: the initial number of foreign cells at 
the time instant 3000. Avg. div t.: the expectation of the waiting 
time between two divisions of the foreign cells. Wins ERS: number 
of wins of the immune system against the foreign cells using the 
ERS model settings. Wins CRS: number of wins of the immune 
system against the foreign cells using the CRS model settings. 


Foreigns 

Avg. div. t. 

Wins ERS 

Wins CRS 

300 

50 

448 

310 

350 

40 

75 

66 

350 

50 

307 

225 

350 

70 

487 

400 

400 

50 

134 

135 


5. Conclusion 

First we described arguments that led us to the ERS theoretical model that em¬ 
phasizes the role of self in creating, maintaining and controlling immune responses 
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Figure 8. Normal immune response against a repeated infection 


Number of the cell reproductions 


200 



■ B cell medium reprod 

■ B cell strong reprod 

■ B cell weak reprod 

■ T cell medium reprod 

■ T cell strong reprod 

■ T cell weak reprod 


Figure 9. Division of weak/intermediate/strong kind of T and B 
cells in the simulation showed in Figure [8] 


to self and nonself. Then we discussed the MiStlmm in silico model that was made 
to investigate some important characteristics of immune development, starting from 
conception and ending some time after birth. Finally, results of some computer ex¬ 
periments were discussed. An important part of the latter was the comparison of 
the CRS and ERS theoretical models. In sum, we think that it is likely that evo¬ 
lution preferred adaptive immune systems whose basic mechanism is closer to the 
ERS model than to a CRS model, because ERS gives better results to overcome a 
critical primary infection. 
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